/* Copyright (c) 2021, National University of Defense Technology. All rights reserved.
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 * 
 *     http://www.apache.org/licenses/LICENSE-2.0
 * 
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

#include <cmath>
#include <iostream>
#include "vec_dot.h"
#include "parvec_axpby.h"
#include "par_lobpcg_eigensolver.h"

#ifdef EIGEN_USE_LAPACKE
#include <lapacke.h>
#endif

namespace YHAMG
{

typedef struct
{
	int size[2];
	int ld;
	double* values;
} ColMajor_Matrix;

static inline ColMajor_Matrix ColMajor(int n, int m, int ld, double* values)
{
	ColMajor_Matrix A;
	A.size[0] = n;
	A.size[1] = m;
	A.ld = ld;
	A.values = values;
	return A;
}

static void SetToZeros(const ColMajor_Matrix& A)
{
	int n = A.size[0];
	int m = A.size[1];
	int ld = A.ld;
	double* Av = A.values;

	for (int j = 0; j < m; ++j)
		for (int i = 0; i < n; ++i)
			Av[i + j * ld] = 0.0;
}

static void SetToIdentity(const ColMajor_Matrix& A)
{
	int n = A.size[0];
	int ld = A.ld;
	double* Av = A.values;

	for (int j = 0; j < n; ++j)
	{
		for (int i = 0; i < n; ++i)
			Av[i + j * ld] = 0.0;
		Av[j + j * ld] = 1.0;
	}
}

static void SetDiagonals(const Vector& D, const ColMajor_Matrix& A)
{
	int n = D.size;
	int ld = A.ld;
	double* Dv = D.values;
	double* Av = A.values;

	for (int j = 0; j < n; ++j)
		Av[j + j * ld] = Dv[j];
}

static int CholeskyFactorize(const ColMajor_Matrix& A)
{
	int n = A.size[0];
	int ld = A.ld;
	double* Av = A.values;

	for (int i = 0; i < n; ++i)
	{
		for (int k = 0; k < i; ++k)
			Av[i + i * ld] -= Av[k + i * ld] * Av[k + i * ld];
		if (Av[i + i * ld] <= 0) return 1;
		Av[i + i * ld] = sqrt(Av[i + i * ld]);
		for (int j = i + 1; j < n; ++j)
		{
			for (int k = 0; k < i; ++k)
				Av[i + j * ld] -= Av[k + i * ld] * Av[k + j * ld];
			Av[i + j * ld] /= Av[i + i * ld];
		}
	}

	return 0;
}

static inline int GeneralizedEigenSolve(const ColMajor_Matrix& A, const ColMajor_Matrix& B, const Vector& D)
{
#ifdef EIGEN_USE_LAPACKE
	return LAPACKE_dsygv(LAPACK_COL_MAJOR, 1, 'V', 'U', A.size[0], A.values, A.ld, B.values, B.ld, D.values);
#endif
	return 1;
}

static void ParMultiVec_Copy(const Par_MultiVector& X, const Par_MultiVector& Y)
{
	int m = X.num_vectors;
	for (int j = 0; j < m; ++j)
		Y(j).Copy(X(j));
}

static void ParMultiVec_CopyI(const Par_MultiVector& X, const Par_MultiVector& Y, const bool* mask)
{
	int m = X.num_vectors;
	for (int j = 0, k = 0; j < m; ++j)
		if (mask[j]) Y(k++).Copy(X(j));
}

static void ParMultiVec_RDU(const Par_MultiVector& X, const ColMajor_Matrix& U)
{
	int m = X.num_vectors;
	int ld = U.ld;
	double* Uv = U.values;
	for (int j = 0; j < m; ++j)
	{
		int i = 0;
		for ( ; i + 1 < j; i += 2)
			X(j).Add2Scaled(-Uv[i + j * ld], X(i), -Uv[i + 1 + j * ld], X(i + 1));
		if (i + 1 == j)
			X(j).AddScaled(-Uv[i + j * ld], X(i));
		X(j).Scale(1.0 / Uv[j + j * ld]);
	}
}

static void ParMultiVec_XAPBY(const Par_MultiVector& X, const ColMajor_Matrix& A, double b, const Par_MultiVector& Y, const Par_MultiVector& W)
{
	int n = X.num_vectors;
	int m = W.num_vectors;
	int ld = A.ld;
	double* Av = A.values;
	for (int j = 0; j < m; ++j)	
	{
		ParVec_AXPBY(b, Y(j), 0.0, W(j), W(j));
		int i = 0;
		for ( ; i + 1 < n; i += 2)
			W(j).Add2Scaled(Av[i + j * ld], X(i), Av[i + 1 + j * ld], X(i + 1));
		if (i + 1 == n)
			W(j).AddScaled(Av[i + j * ld], X(i));
	}
}

static void ParMultiVec_XDPBY(const Par_MultiVector& X, const Vector& D, double b, const Par_MultiVector& Y, const Par_MultiVector& W)
{
	int m = D.size;
	double* Dv = D.values;
	for (int j = 0; j < m; ++j)
		ParVec_AXPBY(Dv[j], X(j), b, Y(j), W(j));
}

static void ParMultiVec_Dot(const Par_MultiVector& X, const Par_MultiVector& Y, const ColMajor_Matrix& A)
{
	int n = X.num_vectors;
	int m = Y.num_vectors;
	int ld = A.ld;
	double* Av = A.values;

	double* temp = new double[n * m];

	for (int i = 0; i < n; ++i)
		for (int j = 0; j < m; ++j)
			temp[i + j * n] = Vec_Dot(X.Loc(i), Y.Loc(j));

	MPI_Allreduce(MPI_IN_PLACE, temp, n * m, MPI_DOUBLE, MPI_SUM, X.comm);

	for (int i = 0; i < n; ++i)
		for (int j = 0; j < m; ++j)
			Av[i + j * ld] = temp[i + j * n];

	delete[] temp;
}

static void ParMultiVec_DotU(const Par_MultiVector& X, const Par_MultiVector& Y, const ColMajor_Matrix& U)
{
	int m = X.num_vectors;
	int ld = U.ld;
	double* Uv = U.values;

	double* temp = new double[(m + 1) * m / 2];

	for (int i = 0; i < m; ++i)
		for (int j = i; j < m; ++j)
			temp[i + (j + 1) * j / 2] = Vec_Dot(X.Loc(i), Y.Loc(j));
	
	MPI_Allreduce(MPI_IN_PLACE, temp, (m + 1) * m / 2, MPI_DOUBLE, MPI_SUM, X.comm);

	for (int i = 0; i < m; ++i)
		for (int j = i; j < m; ++j)
			Uv[i + j * ld] = temp[i + (j + 1) * j / 2];

	delete[] temp;
}

static void ParMultiVec_DotD(const Par_MultiVector& X, const Par_MultiVector& Y, const Vector& D)
{
	int m = X.num_vectors;
	double* Dv = D.values;

	for (int j = 0; j < m; ++j)
		Dv[j] = Vec_Dot(X.Loc(j), Y.Loc(j));

	MPI_Allreduce(MPI_IN_PLACE, Dv, m, MPI_DOUBLE, MPI_SUM, X.comm);
}

static void Vec_Sqrt(const Vector& x)
{
	int n = x.size;
	double* xv = x.values;
	for (int i = 0; i < n; ++i)
		xv[i] = sqrt(xv[i]);
}

static int LOBPCG_CheckResidual(const Vector& res, double tol, bool* mask = 0)
{
	int m = res.size;
	double* resv = res.values;
	int cnt = 0;

	if (mask)
	{
		for (int i = 0; i < m; ++i)
			if (mask[i] = resv[i] > tol) ++cnt;
	}
	else
	{
		for (int i = 0; i < m; ++i)
			if (resv[i] > tol) ++cnt;
	}
	
	return cnt;
}

Par_LOBPCG_EigenSolver::Par_LOBPCG_EigenSolver(const Par_Operator& _A)
	: A(&_A),
	B(0),
	T(0),
	Y(0),
	MaxIters(500),
	Tolerance(1.0e-08),
	PrintStats(0)
{
}

Par_LOBPCG_EigenSolver::Par_LOBPCG_EigenSolver(const Par_Operator& _A, const Par_MultiVector& _Y)
	: A(&_A),
	B(0),
	T(0),
	Y(&_Y),
	MaxIters(500),
	Tolerance(1.0e-08),
	PrintStats(0)
{
}

Par_LOBPCG_EigenSolver::Par_LOBPCG_EigenSolver(const Par_Operator& _A, const Par_Operator& _T)
	: A(&_A),
	B(0),
	T(&_T),
	Y(0),
	MaxIters(500),
	Tolerance(1.0e-08),
	PrintStats(0)
{
}

Par_LOBPCG_EigenSolver::Par_LOBPCG_EigenSolver(const Par_Operator& _A, const Par_Operator& _T, const Par_MultiVector& _Y)
	: A(&_A),
	B(0),
	T(&_T),
	Y(&_Y),
	MaxIters(500),
	Tolerance(1.0e-08),
	PrintStats(0)
{
}

Par_LOBPCG_EigenSolver::Par_LOBPCG_EigenSolver(const Par_Operator& _A, const Par_Operator& _B, const Par_Operator& _T)
	: A(&_A),
	B(&_B),
	T(&_T),
	Y(0),
	MaxIters(500),
	Tolerance(1.0e-08),
	PrintStats(0)
{
}

Par_LOBPCG_EigenSolver::Par_LOBPCG_EigenSolver(const Par_Operator& _A, const Par_Operator& _B, const Par_Operator& _T, const Par_MultiVector& _Y)
	: A(&_A),
	B(&_B),
	T(&_T),
	Y(&_Y),
	MaxIters(500),
	Tolerance(1.0e-08),
	PrintStats(0)
{
}

Par_LOBPCG_EigenSolver::Par_LOBPCG_EigenSolver(const Par_Operator& _A, int max_Iters, double tolerance, int print_stats)
	: A(&_A),
	B(0),
	T(0),
	Y(0),
	MaxIters(max_Iters),
	Tolerance(tolerance),
	PrintStats(print_stats)
{
}

Par_LOBPCG_EigenSolver::Par_LOBPCG_EigenSolver(const Par_Operator& _A, const Par_MultiVector& _Y, int max_Iters, double tolerance, int print_stats)
	: A(&_A),
	B(0),
	T(0),
	Y(&_Y),
	MaxIters(max_Iters),
	Tolerance(tolerance),
	PrintStats(print_stats)
{
}

Par_LOBPCG_EigenSolver::Par_LOBPCG_EigenSolver(const Par_Operator& _A, const Par_Operator& _T, int max_Iters, double tolerance, int print_stats)
	: A(&_A),
	B(0),
	T(&_T),
	Y(0),
	MaxIters(max_Iters),
	Tolerance(tolerance),
	PrintStats(print_stats)
{
}

Par_LOBPCG_EigenSolver::Par_LOBPCG_EigenSolver(const Par_Operator& _A, const Par_Operator& _T, const Par_MultiVector& _Y, int max_Iters, double tolerance, int print_stats)
	: A(&_A),
	B(0),
	T(&_T),
	Y(&_Y),
	MaxIters(max_Iters),
	Tolerance(tolerance),
	PrintStats(print_stats)
{
}

Par_LOBPCG_EigenSolver::Par_LOBPCG_EigenSolver(const Par_Operator& _A, const Par_Operator& _B, const Par_Operator& _T, int max_Iters, double tolerance, int print_stats)
	: A(&_A),
	B(&_B),
	T(&_T),
	Y(0),
	MaxIters(max_Iters),
	Tolerance(tolerance),
	PrintStats(print_stats)
{
}

Par_LOBPCG_EigenSolver::Par_LOBPCG_EigenSolver(const Par_Operator& _A, const Par_Operator& _B, const Par_Operator& _T, const Par_MultiVector& _Y, int max_Iters, double tolerance, int print_stats)
	: A(&_A),
	B(&_B),
	T(&_T),
	Y(&_Y),
	MaxIters(max_Iters),
	Tolerance(tolerance),
	PrintStats(print_stats)
{
}

void Par_LOBPCG_EigenSolver::operator()(const Par_MultiVector& X, const Vector& lambda, int& iter, const Vector& res) const
{
	MPI_Comm comm = A->comm;

	int rank;
	MPI_Comm_rank(comm, &rank);

	int n = A->OutSize();
	int m = X.num_vectors;
	int m1 = m;

	int my = 0;
	if (Y) my = Y->num_vectors;

	Par_MultiVector R(comm), W(comm), P(comm), Y1(comm), BY(comm), AX(comm), BX(comm), AW(comm), BW(comm), AP(comm), BP(comm), TempVec(comm);
	ColMajor_Matrix GramA, GramB, GramY, YBY, YBX, YBW, XAX, XBX, WAW, WBW, PAP, PBP, XAW, XBW, XAP, XBP, WAP, WBP, Yx, Yw, Yp;
	Vector GramLamdba;

	R.Allocate(m, n);
	W.Allocate(m, n);
	P.Allocate(m, n);
	Y1.Allocate(my, n);
	BY.Allocate(my, n);
	AX.Allocate(m, n);
	BX.Allocate(m, n);
	AW.Allocate(m, n);
	BW.Allocate(m, n);
	AP.Allocate(m, n);
	BP.Allocate(m, n);
	TempVec.Allocate(m, n);
	GramLamdba.Resize(3 * m);

	double* GramAv = new double[9 * m * m];
	double* GramBv = new double[9 * m * m];
	double* GramYv = new double[my * (my + 2 * m)];
	
	if (PrintStats && rank == 0)
	{
		std::cout << "Iter\tResidual\n";
		std::cout << "----------------------------------------------------------------\n";
	}

	GramA = ColMajor(3 * m, 3 * m, 3 * m, GramAv);
	GramB = ColMajor(3 * m, 3 * m, 3 * m, GramBv);
	GramY = ColMajor(my, my + 2 * m, my, GramYv);
	
	XAX = ColMajor(m, m, 3 * m, GramAv);
	XBX = ColMajor(m, m, 3 * m, GramBv);
	YBY = ColMajor(my, my, my, GramYv);
	YBX = ColMajor(my, m, my, GramYv + my * my);

	SetToZeros(GramA);
	SetToZeros(GramB);

	bool* mask = new bool[m];
	bool ExplicitGramFlag = 0;

	iter = 0;
	int iter1 = 0;

	if (Y)
	{
		ParMultiVec_Copy(*Y, Y1);

		if (B)
			B->Apply(Y1, BY);
		else
			ParMultiVec_Copy(Y1, BY);

		ParMultiVec_DotU(BY, Y1, YBY);

		if (CholeskyFactorize(YBY))
			goto breakdown;

		ParMultiVec_RDU(Y1, YBY);

		if (B)
			ParMultiVec_RDU(BY, YBY);
		else
			ParMultiVec_Copy(Y1, BY);

		ParMultiVec_Dot(BY, X, YBX);
		ParMultiVec_XAPBY(Y1, YBX, -1.0, X, X);
	}

	if (B)
		B->Apply(X, BX);
	else
		ParMultiVec_Copy(X, BX);

	ParMultiVec_DotU(BX, X, XBX);

	if (CholeskyFactorize(XBX))
		goto breakdown;

	ParMultiVec_RDU(X, XBX);
	
	if (B) 
		ParMultiVec_RDU(BX, XBX);
	else
		ParMultiVec_Copy(X, BX);
	
	SetToIdentity(XBX);

	A->Apply(X, AX);
	ParMultiVec_DotU(AX, X, XAX);

	if (GeneralizedEigenSolve(XAX, XBX, lambda))
		goto breakdown;

	ParMultiVec_XAPBY(X, XAX, 0.0, TempVec, TempVec);
	ParMultiVec_Copy(TempVec, X);
	ParMultiVec_XAPBY(AX, XAX, 0.0, TempVec, TempVec);
	ParMultiVec_Copy(TempVec, AX);

	if (B)
	{
		ParMultiVec_XAPBY(BX, XAX, 0.0, TempVec, TempVec);
		ParMultiVec_Copy(TempVec, BX);
	}
	else
		ParMultiVec_Copy(X, BX);

	while (1)
	{
		R.num_vectors = m;

		ParMultiVec_XDPBY(BX, lambda, -1.0, AX, R);

		ParMultiVec_DotD(R, R, res);
		Vec_Sqrt(res);

		m1 = LOBPCG_CheckResidual(res, Tolerance, mask);

		if (PrintStats && rank == 0)
		{
			std::cout << iter;
			for (int j = 0; j < m; ++j)
			{
				if (mask[j])
					std::cout << "\t" << res.values[j] << "\n";
				else
					std::cout << "\t" << res.values[j] << " *\n";
			}
		}

		if (m1 == 0 || iter == MaxIters) break;

		if (!ExplicitGramFlag)
			ExplicitGramFlag = m1 < m || LOBPCG_CheckResidual(res, 1.0e-08) < m;

		ParMultiVec_CopyI(R, R, mask);
		R.num_vectors = m1;
		W.num_vectors = m1;
		AW.num_vectors = m1;
		BW.num_vectors = m1;

		XAX = ColMajor(m, m, 3 * m, GramAv);
		XAW = ColMajor(m, m1, 3 * m, GramAv + 3 * m * m);
		WAW = ColMajor(m1, m1, 3 * m, GramAv + m + 3 * m * m);

		XBX = ColMajor(m, m, 3 * m, GramBv);
		XBW = ColMajor(m, m1, 3 * m, GramBv + 3 * m * m);
		WBW = ColMajor(m1, m1, 3 * m, GramBv + m + 3 * m * m);

		if (T)
			T->Apply(R, W);
		else
			ParMultiVec_Copy(R, W);

		YBW = ColMajor(my, m1, my, GramYv + my * (my + m));

		if (Y)
		{
			ParMultiVec_Dot(BY, W, YBW);
			ParMultiVec_XAPBY(Y1, YBW, -1.0, W, W);
		}

		ParMultiVec_Dot(BX, W, XBW);
		ParMultiVec_XAPBY(X, XBW, -1.0, W, W);

		if (B)
			B->Apply(W, BW);
		else
			ParMultiVec_Copy(W, BW);
			
		ParMultiVec_DotU(BW, W, WBW);

		if (CholeskyFactorize(WBW))
			goto breakdown;

		ParMultiVec_RDU(W, WBW);

		if (B)
			ParMultiVec_RDU(BW, WBW);
		else
			ParMultiVec_Copy(W, BW);

		A->Apply(W, AW);
		ParMultiVec_Dot(AX, W, XAW);
		ParMultiVec_DotU(AW, W, WAW);

		if (ExplicitGramFlag)
		{
			ParMultiVec_DotU(AX, X, XAX);
			ParMultiVec_DotU(BX, X, XBX);
			ParMultiVec_DotU(BW, W, WBW);
			ParMultiVec_Dot(BX, W, XBW);
		}
		else
		{
			SetToZeros(XAX);
			SetDiagonals(lambda, XAX);
			SetToIdentity(XBX);
			SetToIdentity(WBW);
			SetToZeros(XBW);
		}

restart:
		if (iter1 > 0)
		{
			XAP = ColMajor(m, m1, 3 * m, GramAv + 3 * m * (m + m1));
			WAP = ColMajor(m1, m1, 3 * m, GramAv + m + 3 * m * (m + m1));
			PAP = ColMajor(m1, m1, 3 * m, GramAv + m + m1 + 3 * m * (m + m1));
			
			XBP = ColMajor(m, m1, 3 * m, GramBv + 3 * m * (m + m1));
			WBP = ColMajor(m1, m1, 3 * m, GramBv + m + 3 * m * (m + m1));
			PBP = ColMajor(m1, m1, 3 * m, GramBv + m + m1 + 3 * m * (m + m1));
			
			GramA = ColMajor(m + m1 + m1, m + m1 + m1, 3 * m, GramAv);
			GramB = ColMajor(m + m1 + m1, m + m1 + m1, 3 * m, GramBv);
			GramLamdba.size = m + m1 + m1;
			
			ParMultiVec_CopyI(P, P, mask);
			ParMultiVec_CopyI(AP, AP, mask);
			ParMultiVec_CopyI(BP, BP, mask);
			P.num_vectors = m1;
			AP.num_vectors = m1;
			BP.num_vectors = m1;

			ParMultiVec_DotU(BP, P, PBP);

			if (CholeskyFactorize(PBP))
			{
				iter1 = 0;
				goto restart;
			}

			ParMultiVec_RDU(P, PBP);
			ParMultiVec_RDU(AP, PBP);

			if (B)
				ParMultiVec_RDU(BP, PBP);
			else
				ParMultiVec_Copy(P, BP);
			
			ParMultiVec_Dot(AX, P, XAP);
			ParMultiVec_Dot(AW, P, WAP);
			ParMultiVec_DotU(AP, P, PAP);
			
			ParMultiVec_Dot(BX, P, XBP);
			ParMultiVec_Dot(BW, P, WBP);

			if (ExplicitGramFlag)
				ParMultiVec_DotU(BP, P, PBP);
			else
				SetToIdentity(PBP);

			if (GeneralizedEigenSolve(GramA, GramB, GramLamdba))
			{
				iter1 = 0;
				goto restart;
			}

			GramLamdba.size = m;
			Yx = ColMajor(m, m, 3 * m, GramAv);
			Yw = ColMajor(m1, m, 3 * m, GramAv + m);
			Yp = ColMajor(m1, m, 3 * m, GramAv + m + m1);

			lambda.Copy(GramLamdba);

			ParMultiVec_XAPBY(P, Yp, 0.0, TempVec, TempVec);
			ParMultiVec_XAPBY(W, Yw, 1.0, TempVec, TempVec);
			P.num_vectors = m;
			ParMultiVec_Copy(TempVec, P);
			
			ParMultiVec_XAPBY(X, Yx, 1.0, TempVec, TempVec);
			ParMultiVec_Copy(TempVec, X);
			
			ParMultiVec_XAPBY(AP, Yp, 0.0, TempVec, TempVec);
			ParMultiVec_XAPBY(AW, Yw, 1.0, TempVec, TempVec);
			AP.num_vectors = m;
			ParMultiVec_Copy(TempVec, AP);

			ParMultiVec_XAPBY(AX, Yx, 1.0, TempVec, TempVec);
			ParMultiVec_Copy(TempVec, AX);

			if (B)
			{
				ParMultiVec_XAPBY(BP, Yp, 0.0, TempVec, TempVec);
				ParMultiVec_XAPBY(BW, Yw, 1.0, TempVec, TempVec);
				BP.num_vectors = m;
				ParMultiVec_Copy(TempVec, BP);

				ParMultiVec_XAPBY(BX, Yx, 1.0, TempVec, TempVec);
				ParMultiVec_Copy(TempVec, BX);
			}
			else
			{
				ParMultiVec_Copy(P, BP);
				ParMultiVec_Copy(X, BX);
			}
		}
		else
		{
			GramA = ColMajor(m + m1, m + m1, 3 * m, GramAv);
			GramB = ColMajor(m + m1, m + m1, 3 * m, GramBv);
			GramLamdba.size = m + m1;

			if (GeneralizedEigenSolve(GramA, GramB, GramLamdba))
				goto breakdown;

			GramLamdba.size = m;
			Yx = ColMajor(m, m, 3 * m, GramAv);
			Yw = ColMajor(m1, m, 3 * m, GramAv + m);

			lambda.Copy(GramLamdba);

			ParMultiVec_XAPBY(W, Yw, 0.0, TempVec, TempVec);
			P.num_vectors = m;
			ParMultiVec_Copy(TempVec, P);
			
			ParMultiVec_XAPBY(X, Yx, 1.0, TempVec, TempVec);
			ParMultiVec_Copy(TempVec, X);

			ParMultiVec_XAPBY(AW, Yw, 0.0, TempVec, TempVec);
			AP.num_vectors = m;
			ParMultiVec_Copy(TempVec, AP);

			ParMultiVec_XAPBY(AX, Yx, 1.0, TempVec, TempVec);
			ParMultiVec_Copy(TempVec, AX);

			if (B)
			{
				ParMultiVec_XAPBY(BW, Yw, 0.0, TempVec, TempVec);
				BP.num_vectors = m;
				ParMultiVec_Copy(TempVec, BP);

				ParMultiVec_XAPBY(BX, Yx, 1.0, TempVec, TempVec);
				ParMultiVec_Copy(TempVec, BX);
			}
			else
			{
				ParMultiVec_Copy(P, BP);
				ParMultiVec_Copy(X, BX);
			}
		}
		
		++iter;
		++iter1;
	}

breakdown:
	if (PrintStats && rank == 0)
	{
		std::cout << "----------------------------------------------------------------\n";
		std::cout << "Iterations: " << iter << "\n";
		std::cout << "Converged Eigenpairs: " << m - m1 << "/" << m << "\n";
	}

	delete[] GramAv;
	delete[] GramBv;
	delete[] GramYv;
	delete[] mask;
}

}
